Homogenization of nonlocal wire metamaterial via a renormalization approach 
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It is well known that defining a local refractive index for a metamaterial requires that the wave- 
length be large with respect to the scale of its microscopic structure (generally the period). However, 
the converse does not hold. There are simple structures, such as the infinite, perfectly conducting 
wire medium, which remain non-local for arbitrarily large wavelength-to-period ratios. In this work 
we extend these results to the more realistic and relevant case of finite wire media with finite 
conductivity. In the quasi-static regime the metamaterial is described by a non-local permittivity 
which is obtained analytically using a two-scale renormalization approach. Its accuracy is tested 
and confirmed numerically via full vector 3D finite element calculations. Moreover, finite wire media 
exhibit large absorption with small reflection, while their low fill factor allows considerable freedom 
to control other characteristics of the metamaterial such as its mechanical, thermal or chemical 
robustness. 
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The effective medium theory of artificial metallo- 
dielectric structures goes back to the beginning of the 
20th century, with the work of Maxwell-Garnetti and 
Wiener-. These, and subsequent effective medium theo- 
ries focused on disordered media where only partial in- 
formation on the microscopic structure was available. A 
major step forward was made with the work of Kock, in 
the 1940s 3 -. This time Lorentz theory 4 ^ was used to de- 
sign artificial effective media, in a bottom up fashion, as 
an array of scatterers. In the 1970s more mathematically 
sophisticated methods emerged, where instead of seeking 
a limiting effective medium (equivalent in some suitably 
defined sense to the structure of interest), one obtains a 
limiting equation system^, for the macroscopic electro- 
magnetic field in a given structure^—. 

In recent years, the advent of negative index metama- 
terials and composites has led to increased interest in 
effective medium theories. The most popular by far is of 
course the Lorentz theory approach, it being the most 
accessible and intuitively appealing^ 2 -. However, the use- 
fulness of Lorentz theory is much diminished when one is 
interested in materials where the size of objects is much 
larger than the distances separating them, or materials 
which are strongly non-local, or in which the scatterers 
are strongly coupled, leading to behavior of a collective 
natur e 13 ' 14 . Contrary to common intuition, non-local be- 
havior persists, in certain structures, even when the wa- 
velength is much larger than the characteristic scale of 
the structure ; an excellent example is the wire medium 
studied by Belov et alJ^— . In these situations the Lo- 
rentz model is no longer useful and more sophisticated 
techniques are required. 

In this work we test and illustrate, for the first time, 
an effective medium model of the finite conductivity fi- 
nite wire medium (the "bed-of-nails" structure, Fig. QJ 
based on a two-scale renormalization approach. Instead 
of letting the wavelength tend to infinity, as customary in 
effective medium theories, we keep it fixed, and let other 
geometrical parameters tend to zero. The advantage of 



this approach is that it leaves us the possibility of kee- 
ping some of the geometrical parameters fixed (in this 
case the wire length L), leading to a new type of partial 
homogenization scheme. To put it less formally, we would 
like to homogenize while keeping the thickness fixed with 
respect to the wavelength, which prevents us from let- 
ting A tend to infinity, so the only remaining option is to 
make all the other dimensions (the wire radius r and the 
period d) tend to zero. 

Unlike common practice in much of the metamaterials 
literature, we include a detailed discussion of the mo- 
del's domain of applicability, so that an engineer may be 
able to quickly and efficiently decide whether this kind 
of structure may be useful for a given purpose. 

The structure under study is a square biperiodic array 
of thin wires, of length L, radius r and conductivity a. 
We note the period d and the wavelength A. The renor- 
malization (depicted in Fig. [T]) involves a limiting process 
whereby the three quantities : r, d and 1/er tend simul- 
taneously to zero. The parameter governing the limiting 
process is noted r\ = d, the period. The asymptotics of 
the other two parameters, a and r, with respect to 7/ are 
described by fixed parameters k and 7 according to the 
following relations : 
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where u> is the angular frequency of the electromagne- 
tic field. In other words the conductivity is renormalized 

inversely to the fill factor 6 V = -^r, while the radius is 
renormalized such that the expression r] 2 log(^-) remains 
constant. 

While these expressions may at first seem obscure, they 
have simple intuitive interpretations. The first requires 
the current density to remain constant during the renor- 
malization. Notice that n is nothing other than the vo- 



2 



Renormalization 




Figure 1: The bed-of-nails structure and the renormalization process. The conducting fibers occupy a region Q £ R 2 , are 
oriented in the z direction, and the structure is periodic in the xy plane. Two renormalized structures are shown, corresponding 
to 771 and 772 respectively, with 771 > 7/2, d m > d V2 , a vi < o m and r vi /d vi > r^/d^ (see Eqs. \T\ and [5J. The real physical 
structure corresponds by definition to 77 = 1 : d v= i — d. The length L and the wavelength A remain fixed, i.e., we are 
homogenizing in the xy plane only. 



lume average of the imaginary part of the permittivity. 
Also, recall that the static admittance per unit length of 
a circular wire is given by 
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and that the number of wires per unit area is given by 
I/77 2 . The second expression requires the average inter- 
nal capacitance of the wires to remain constant during 
renormalization. This feature is known to be essential for 
their asymptotic behavior (see, for instance Refs ; 18 ' 19 ). 
One may object that the expression on the right side of 
Eq.[2]is valid for infinitely long wires, whereas we are wor- 
king with wires of finite length. Indeed, as shown below, 
the model fails for short wires (comparable to the period) , 
but that configuration is best treated with the Lorentz 
approach anyway^, placing it outside our present scope. 

The essential quantities in the rescaling process are 
therefore the geometric quantities r v , 77, the material 
quantity a v and the field quantities E v and H n . To these 
one must also add a quantity characterizing the all im- 
portant electric field in the wires. This is noted F v , it is 
non-zero only inside the wires, and is given by 
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F v has the units of electric field, and in the microscopic, 
inhomogeneous picture it is clearly proportional to the 
current density. In the macroscopic, homogeneous pic- 
ture, however, it will correspond to the polarization den- 
sity P. More precisely we have lim^o F v — P/ieo. 



The question to be answered now becomes : what hap- 
pens in the limit 77 — 5- ? The answer is that the fields 
converge (in a precise sense described in Ref<2i) to the 
unique solution of the following system : 
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dz 



= iuj/ioH 
= —iojEo(E + j^z) 
= -2Trje E z , ze [-L/2,L/2] 
= 0, z e {-L/2.L/2} 

(4) 

Before solving the system, let us first see what it tells us 
on a more intuitive level. 

All field quantities above are effective, homogeneous 
quantities, which have meaning when the wires have been 
replaced with a homogeneous effective medium with an 
electric polarization density equal to P. The equation 
which gives P is an inhomogeneous Helmholtz equation 
where the source term is given by the z component of 
the electric field E z . The polarization satisfies Neumann 
conditions at the upper and lower interfaces of the slab. 
It is not in general continuous there because Maxwell's 
equations impose the continuity of the normal component 
of the displacement field D = EqE + P ; consequently, 
any jump in E must be canceled by an equivalent jump 
in P/eq. The dependence of P on E, i.e., the constitu- 
tive relation, takes the form of an integral. In this case we 
are dealing with a one-dimensional inhomogeneous Helm- 
holtz equation, but this situation is slightly complicated 
by the fact that it is valid on a bounded domain only 
(the thickness L of the slab) . The polarization field has 
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the form 
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g(z,z )E z (x,z)dz (5) 



where g(z, z ) is the Green function of the Helmholtz 
operator on the bounded domain z 6 (— -§,-§■)• It takes 
the form (see Appendix A) 
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where K 2 = feg + ^p, z < = mw.(z,Zo) and z > — 
max(z,zo)- Relation [5] is clearly a non-local constitutive 
relation because the value of the polarization field at a 
position zq depends on values of the electric field at po- 
sitions different from zq. 

When the imaginary part of K is large the integral 
above drops off quickly. In the limit of small conductivity 
(and hence small k), the polarization becomes local for 
sufficiently large wavelengths. In the opposite limit, for 
infinite conductivity and infinitely long wires the integral 
covers all space (in the z direction) and the material is 
non-local, even in the long- wavelength regime. In fact this 
can be seen immediately by doing a Fourier transform on 
the third equation of system (with k — » oo) : 
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This is consistent with the findings of Belov et al»i£~— . 

Until now, the discussion has been independent of the 
actual shape of the domain (Fig. [TJ. From this point 
on, however, for purposes of illustration we specialize to 
the case Vt = R 2 , which is an infinite two dimensional 
bed-of-nails, of thickness L, period d, wire radius r and 
conductivity a. The effective medium is therefore a ho- 
mogeneous slab parallel to the xy plane and of thickness 
L. 
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Figure 2: Transmission (solid), reflection (dot-dashed) and 
absorption (dashed) efficiency curves comparing the finite ele- 
ment solution (dot markers) and the effective medium solution 
(no markers) as a function of angle of incidence. The struc- 
ture has a conductivity a — 8(f2m) _1 , period d — 0.01m, and 
dimensionless parameters L/d = 120, \/d — 20, r/d = 0.1, 
and S/d = 4.6. Computational constraints forced us to use 
a very coarse mesh, which explains the approximate nature 
of the energy conservation (x markers) of the finite element 
model. 
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Figure 3: Square of the current density for the effective me- 
dium solution of Eqs. [9] (dashed) and the finite element solu- 
tion (solid) as a function of position within the bed-of-nails 
structure (which is positioned in z G (0, L)). The structure 
is the same as in Fig. [2] illuminated at an angle of incidence 
9 = 40° from the top. 



Numerical results 

We now proceed to test the homogeneous model by 
comparing it with 3D full vector simulations of the struc- 
ture, i.e. we must compare the reflection, transmission 
and absorption coefficients and the current distribution 
of the homogeneous problem with those of the original 
bed-of-nails metamaterial. The solution to the homoge- 
neous problem is obtained by integrating system 0] as 
described in Appendix B. 

The 3D full vector simulations of the bed-of-nails meta- 
material were done using the Comsol Multiphysics finite 



element method 22 software package. The periodicity was 
implemented using Floquet-Bloch conditions 2 ^ in the two 
periodic directions (x and y), and absorbing Perfectly 
Matched Layers 2 ^ in the positive and negative z direc- 
tions. The linearity of the materials in the structure was 
used to treat the incident field as a localised source within 
the obstacle, as detailed in Ref i 25 ' 26 . The Comsol/Matlab 
scripts of the models used to produce the figures below 
are available as online support material for the readers' 
convenience. 

Figures [5] and [3] show good agreement between the ef- 
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fective medium model and the finite element simulation. 
Note that the current density behavior near the bounda- 
ries differs between the effective medium model and the 
finite element model. This is due to the fact that in the 
macroscopic, homogeneous scenario, one speaks of a po- 
larization field obeying Neumann boundary conditions, 
as discussed above. In the microscopic scenario however, 
we have a free conductor carrying current induced by 
an external electric field. Since in our geometry at the 
given wavelength the capacitance of the wire endpoints 
is very small, the accumulation of charge will be corres- 
pondingly small, leading to an almost continuous normal 
component of the electric field (and therefore also cur- 
rent). Numerically, it seems as if the current goes to zero 
at the wire endpoints, even though this is not strictly 
exact. Nevertheless, since in the homogeneous limit the 
boundary condition of the current is of Neumann type, 
the convergence of the renormalization process is clearly 
non-uniform near the boundaries. This provides an addi- 
tional explanation for requiring long wires ; we want the 
effect of the boundaries to be small. 

It must also be pointed out that the parameters of the 
particular structure chosen for the illustration in Figs. [5] 
and[3]were forced upon us by practical constraints : finite 
element meshing of thin long circular wires requires very 
large amounts of computer memory and time. Simulation 
of wires thinner than r/d = 0.05 is prohibitive. Conse- 
quently, in order to explore a wider domain of the pa- 
rameter space, we have taken advantage of the fact that 
the structures we are interested in have r <C d and 8 3> r. 
Such thin conducting structures can be simulated much 
more efficiently as lines of zero thickness^ (i.e. edges, in 
the finite element formulation) carying current and ex- 
hibiting an equivalent linear impedance. This approach 
gives excellent results with a fraction of the computing 
power, and enables us to model realistic structures that 
would otherwise be inaccessible. 

For instance, Figs. Q] and [S] show the results of calcu- 
lations for a structure of Toray T300®carbon fibers^ 
with a conductivity of : a = 5.89 • 10 4 (f2m) _1 and a 
radius of 3.5 microns. The wires have an aspect ratio 
L/r = 2.28 x 10 5 , which is far beyond what would have 
been accessible by meshing the interior of the wires. The 
finite element model of Fig. [5] (curves with markers) , in 
which the interior of the wires is meshed, is a problem 
with approx. 2.8 million degrees of freedom, which re- 
quires at least 42 Gigabytes of available RAM to solve. 
By comparison, the model of Fig. U (curves with mar- 
kers) , in which the wires are modeled as current carrying 
edges, is a problem of approx. 62 thousand degrees of 
freedom, which requires less than one Gigabyte of avai- 
lable RAM and can therefore be solved on any sufficiently 
recent desktop computer. 

Figures [21 El HI and 03 illustrate the behavior which is 
typical of the model. The agreement remains good up to 
high incidence angles, and over a large wavelength do- 
main (Fig. [7]). The structure is transparent in normal 
incidence. For increasingly oblique angles of incidence 
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Figure 4: Transmission (solid), reflection (dot-dashed) and 
absorption (dashed) efficiency curves comparing the finite 
element solution (dot markers) and the effective medium 
(no markers) as a function of angle of incidence. The wire 
conductivity is that of Toray T300® carbon fibers a = 
5.89-10 4 (fim) -1 . The structure has period d = 0.01m, and di- 
mensionless parameters L/d = 80, X/d = 20, r/d = 3.5 TO -4 , 
and S/r — 15. Energy conservation of the finite element mo- 
del ( x markers) is respected to within better than one percent 
for most angles of incidence. The departure around 80° is ex- 
plained by the poor performance of the PML absorbing layers 
when close to grazing incidence. 
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Figure 5: Square of the current density for the effective me- 
dium model (dashed) and the finite element solution (solid) 
as a function of position within the slab (which is positioned 
in z E (0, L)). The structure is the same as in Fig. [31 illumina- 
ted at an angle of incidence 9 = 40° from the top. Note that 
the surface areas under the two curves (in this figure as well 
as Fig. [3} are the same because they are proportional to the 
Joule dissipation rates, which are seen to be equal from Fig. 
[4] (and Fig. [2]) at the given angle of incidence. 



the absorption increases more or less gradually, depen- 
ding on the thickness L. The reflection is generally low, 
though it increases when approaching grazing incidence. 
The low reflection may be explained by the small radii of 
the wires : their extremities have low capacitance, hence 
they exhibit very little charge accumulation, leading to 
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Figure 6: Transmission (solid), reflection (dot-dashed) and 
absorption (dashed) efficiency curves comparing the finite ele- 
ment solution (dot markers) and the effective medium (no 
markers) as a function of angle of incidence. The structure 
has a conductivity a = lOOO(fim) -1 , period d — 0.01m, and 
dimensionless parameters L/d = 50, \/d = 8, r/d = 0.002, 
and 5/d = 13. The reflection remains low for angles of inci- 
dence of up to 80° even as the Joule absorption reaches almost 
100% for 9 > 60°. Energy conservation is indicated by the x 
markers. 

an almost continuous normal component of the electric 
field. Certain configurations exhibit very low reflection 
for almost all angles of incidence, see Fig. [5] and Fig. [7] 
around A = 1.2m. The current density decreases roughly 
exponentially within the structure due to absorption. 

Domain of validity 

The boundaries of the domain of validity of the model 
are given by four dimensionless parameters : the ratio of 
the skin depth to the radius in the wires S/r, the ratio 
of the wire length to the period L/d, the ratio of the 
wavelength to the period X/d and the ratio of the wire 
radius to the period r/d. 

The skin depth must be larger than the radius, due to 
the fact that the impedance used in defining k (Eq. [3]) 
is the static impedance which differs from the quasistatic 
value by an imaginary inductive term iw/x/87r (see, for 
instance, Ref^). Requiring this term to be negligible is 
equivalent to requiring that S 2 /r 2 1. Moreover, in the 
rescaling process the skindepth/radius ratio is given by 

s JL = x rr 

r n rj V 2ttk 

Since ?y approaches zero in the rescaling process, it is 
natural to expect the homogeneous model to be valid 
when the skindepth is large compared to the radius. 

In addition, recall that the definition of 7 in Eq. [2] 
fixes the capacitance of the wires to the value for thin, 
long wires. Consequently, we expect the model to hold 
for large L/d and for small r/d. To these, we must add 



Figure 7: Transmission (solid), reflection (dot-dashed) and 
absorption (dashed) efficiency curves comparing the finite ele- 
ment solution (dot markers) and the effective medium (no 
markers) as a function of wavelength. Energy conservation 
for the finite element model is labeled with x markers. The 
structure has a conductivity a = 3000(fim) _1 (in the semi- 
conductor domain), period d = 0.01m, and dimensionless pa- 
rameters L/d = 60, r/d = 0.003, and the angle of incidence 
is 8 = 70°. 5/r runs approximately from 4 to 25 from left 
to right over the domain of the plot. The model fails around 
A < 0.1m = lOd. 



the general requirement for all effective medium models : 
the wavelength must be large compared to the period. 

Due to the large (four dimensional) parameter space, 
an exhaustive numerical exploration of the bed-of-nails 
structure is not feasible in a reasonable timeframe. Still, 
our study has made it possible to broadly determine the 
boundaries of the domain of applicability of the effective 
medium model. Roughly, one must have X/d 7 — 12, 
5/r > 4 - 8, L/d > 20 - 30, r/d < 10. Our (a for- 
tiori limited) numerical exploration of the parameter 
space suggests that the skindepth-to-radius ratio is often 
the main limiting factor, particularly when considering 
highly conducting wires. 



Conclusion 

We have tested numerically the effective medium 
theory of the bed-of-nails structure, whose rigorous ma- 
thematical foundation is described in RefiSl. We have 
found good agreement between the transmission, reflec- 
tion and absorption efficiencies between the effective me- 
dium model and a 3D finite element model, for a broad 
range of angles of incidence and wavelengths. The current 
density in the real structure corresponds to the polariza- 
tion current density of the effective medium model. The 
medium is nonlocal, meaning that the polarization field 
depends on the electric field over a region of finite size. 
That dependence is given by Eq. [SJ This nonlocal be- 
havior also means that the permittivity depends on the 
wavevector, so it can no longer be seen, strictly, as a pro- 
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perty of the medium, but rather, as a property of a given 
wave propagating in the structur o 30 i 31 . 

The bed-of-nails structure is a medium exhibiting high 
absorption with low reflection. It requires a very low 
filling fraction of conducting material, but exhibits near 
perfect absorption over a wide range of angles of inci- 
dence, for sufficiently large thicknesses. The low filling 
fraction is useful because it allows the engineer to fill the 
space between the wires with materials satisfying other 
design constraints, such as mass density, or mechanical, 
chemical or thermal robustness. The geometries studied 
here are transparent at normal incidence, but this aspect 
can easily be rectified by slanting the wires by about 20° 
with respect to the upper and lower boundaries. This 
design may therefore be used to obtain a near-perfect 
electromagnetic absorber for all angles of incidence in a 
very straightforward way and with considerable freedom 
in the resulting mechanical, thermal or chemical proper- 
ties of the structure. We are currently exploring more 
elaborate structures which may be modeled by the same 
scaling technique : structures with thin wires in the x 
and/or y directions as well as the z direction, or with 
wires curved helically, leading to a non-trivial magnetic 
constitutive relation in addition to the electric one. 



Appendix A 

We require the Green function for the problem (see 
chapter II of Ref .22) 



p" + c?p = (3E Z 
with 

p'(-L/2)=p(L/2)=0. 



(G) 



For the purpose of this Appendix it is convenient to consi- 
der the structure is positioned between— L/2 and L/2. 
The Green function satisfies the equation 



g" + a 2 g = 5 Zo , z £ [ 
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and may be written : 

g(z,z Q ) = Cu 1 (z < )u 2 (z > ) 
with 

z < = min(z, Zq) 
z > = max(z, zq) 

such that 



(7) 

(8) 



when z £ | — 



L 



Cui(z)u 2 {z ) 



and 



when z £ ^z , — J , g = Cu\(z )u 2 {z) 

Replacing form [8] into Eq. [7] one obtains that g must 
be continuous at Zo, its derivative must have a jump dis- 
continuity of 1, and the two functions u\ and u 2 must be 



sinusoidal of wave constant a : 

ut{z) = Acos(a(z + L/2)) 
u 2 {z) = Bcos(a{z - L/2)). 

By imposing the boundary conditions Eq. |B] we obtain 

u'ii-L/2) = 

u 2 (L/2) = 

and by requiring a jump discontinuity of 1 at Zq we obtain 

1 



ABC 



giving finally : 

g{z,z ) 



i 



a sin(aL) 



a sin(aL) 
cos(a(z< + L/2)) cos(a(z> - L/2)) 
Appendix B 



We now proceed to solve the homogeneous limit sys- 
tem |4] For convenience we position it in z £ (0, L). Since 
we are dealing with a system with translational inva- 
riance, a slab, we can split the problem into two inde- 
pendent polarization cases : TE, where the electric field 
is in the xy plane, and TM, where the magnetic field is 
in the xy plane. However, since we are considering thin 
wires (small volume fraction) the structure will be trans- 
parent to TE waves. We therefore only have to consider 
TM waves. We choose a coordinate system so that the 
plane of incidence is the xz plane, with angle of inci- 
dence 9, in which case our unknowns will be H y and P z . 
The translation invariance allows us to seek solutions of 
the form : 

H y = u(z)e lax 
P z = p(z)e lax 

with : a = kosmO. Inserting these into system @] we ob- 
tain a system of equations for u and p : 

!u"{z) + (fcg — a 2 ) u(z) = autp(z) 

V'\z) + (k 2 + H^a _ 2n 1 )p(z) = 2 -^u(z), z £ [0,L] 

(9) 

with the important boundary conditions : p' = at 
z = and z = L, and u and v! continuous everywhere. 

The objective is now to obtain the transfer matrix T 
of the slab, which relates the field u and its derivative u' 
at the bottom and the top of the slab : 

u(0) 
u'(0) 

Once T is known the reflection and transmission coeffi- 
cients r and t can be obtained immediately from 

2e~^ L , . 

(11) 



u(L) 
u'(L) 



= T 



(10) 



A-B 
A = Tii - i/3T 12 and B 



A-B 
_ T 2 i - i/3T 22 
i(3 



7 



where f3 = ko cos 9 = \/k 2 — a 2 . 

We begin by integrating system [9] Noting S 2 
"^HJ- — 2nj for readability, we rewrite the system as 



;-2 

K 



W"(z) = -MW(z) 



(12) 



where 



and 



W(z) 



u(z) 
p(z) 



M 



P 2 



-acu 

s 2 



The matrix M can be diagonalized M = QDQ 1 with 
D = di&g(K 2 ,K 2 ) so the system IT21 can be rewritten 
Q- l W"(z) = -DQ- x W(z). Since Q is constant and 
known, this can be integrated directly, and the general 
solution is then obtained as a sum of plane waves : 



^Wiz) 



A+ exp(iK u z) 
A+ exp(iK p z) 



A u exp(-iK u z) , . . 
A- exp(-iX p z) l [i6 > 



Once the integration performed, obtaining T is now only 
a matter of algebraic manipulation, u and p are now ex- 
pressed in terms of the elements of the matrix Q and the 
coefficients A+, A~, A+ and A~ . However, recall that 
we are not interested directly in these coefficients, but 
in the matrix T . Since that matrix does not depend di- 
rectly on p the first step is to eliminate the A p s from the 
equation system. This is done by making use of the boun- 
dary conditions. By differentiating the bottom equation 
of system HHl we can obtain p' as 



p = iK u Q 2 i(A+e 



— A,, 



iK u z 



+iK p Q 22 (A+e lK " z - A 



- -%K 
P e 



) 

pZ ). 



Setting this to zero at z = 0, L we can obtain the A p s in 
terms of the A u s. Noting vectors 



An. 



A, t 



we introduce the matrix 

K U Q 21 1 



C = - 



K p Q 22 2i sin(KpL) 

iK u L p — iK p L 



iK u L 



— iKpL 
AKpL 



e -iK u L 
n -iK u L 



so that 



Ap — CA U . 



We are now in a position to express W(z) in terms of A u 
alone. Eq. [T3] can be rewritten 



W{z) = QE(z)Au 
where E(z) is defined as 



(14) 



i.K„ 



E ^ { Cne lK " z + C 21 e- lK " z C 12 e iK " z + C 22 e 



Eq. Q3] contains (within its first row) the expression for 
u. But to obtain the transfer matrix T we also require u' . 
We simply differentiate Eq. [TJ] to obtain 



W'(z) = QE'{z)Au. 



(15) 



By combining the first rows of Eqs. [Ml and [TBI we are 
in a position to construct the matrix G(z) such that 



u(z) 
u'(z) 



G(z)A u . 



By writing this equation at z = and z — L we obtain 



u , {z)J -G(L)G(0) , /(()) 



Comparing with Eq. [TU] we obtain the result we seek, 

T = G(L)G(0)- 1 , 

leading to the reflection and transmission coefficients via 
Eqs. [TTJ The Matlab script of the above manipulations 
is available as online support material for the readers' 
convenience. 

To summarize, we are now capable of modeling a struc- 
ture with a given d, r, a, L at a given incident field wa- 
velength A in the following way. We first obtain the two 
rescaling parameters k and 7 for the given structure using 
Eqs. [TJ and [2] Then, we integrate system [9] to obtain the 
reflection and transmission coefficients. 
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